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Abstract 

Background: Given the complex mechanisms underlying biochemical processes systems biology researchers tend 
to build ever increasing computational models. However, dealing with complex systems entails a variety of 
problems, e.g. difficult intuitive understanding, variety of time scales or non-identifiable parameters. Therefore, 
methods are needed that, at least semi-automatically, help to elucidate how the complexity of a model can be 
reduced such that important behavior is maintained and the predictive capacity of the model is increased. The 
results should be easily accessible and interpretable. In the best case such methods may also provide insight into 
fundamental biochemical mechanisms. 

Results: We have developed a strategy based on the Computational Singular Perturbation (CSP) method which 
can be used to perform a "biochemically-driven" model reduction of even large and complex kinetic ODE systems. 
We provide an implementation of the original CSP algorithm in COPASI (a COmplex PAthway Simulator) and 
applied the strategy to two example models of different degree of complexity - a simple one-enzyme system and 
a full-scale model of yeast glycolysis. 

Conclusion: The results show the usefulness of the method for model simplification purposes as well as for 
analyzing fundamental biochemical mechanisms. COPASI is freely available at http://www.copasi.org. 



1 Background 

Biochemical systems are inherently high dimensional 
due to the large number of interrelated cellular compo- 
nents and processes, the temporal organization of which 
spans time scales of several orders of magnitude. Aiming 
at a comprehensive understanding of the dynamic beha- 
vior of such systems has led to the development of an 
ever increasing number of computational models which 
are in the majority of cases formulated on the basis of 
ordinary differential equations (ODEs) [1]. Even though 
the purpose of computational models is to facilitate 
understanding and analysis of the underlying biochem- 
ical mechanisms, this again becomes cumbersome with 
the growing complexity of modern models. Therefore, it 
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is important to identify those parts of the biochemical 
systems and of the model that are responsible for the 
observed physiological behavior. This necessitates the 
development of methods for the rational simplification 
of computational models and to make them comfortably 
accessible to the community. 

Numerous methods have been developed to simplify 
(bio)chemical reaction networks (see review [2]). For 
biochemical systems many of the reduction methods 
aim at analyzing the steady state behavior either heuris- 
tically [3] or employing mathematical analysis (e.g. sen- 
sitivities [4,5]). Since biochemical systems usually do not 
reside in a steady state time-dependent approaches have 
recently been proposed (see for example [6,7]). Most of 
these use a mathematical analysis of the different time- 
scales occurring in the biochemical systems, e.g. the 
Intrinsic Low-Dimensional Manifolds (ILDM) method 
[8-11] and the Computational Singular Perturbation 
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(CSP) method [12-14]. Apart from the advantage of a 
time-resolved analysis, these methods can provide useful 
insights, such as the support of the detection of fast 
reactions and species as well as the identification of 
potential rate controlling reactions. However, a disad- 
vantage of the above methods is that the reduced mod- 
els are systems of mathematically transformed 
differential or differential algebraic equations (DAE) 
which may not relate one-to-one to biochemical species 
and reactions hampering the biochemical interpretation. 
In contrast, the methods based on steady-state or partial 
equilibrium approximations keep the one-to-one rela- 
tion and are therefore simple to biochemically interpret. 

In this paper, we focus on deriving simplified bio- 
chemical models by discarding fast reactions and spe- 
cies. For this purpose we present a reduction strategy 
which is based on the CSP algorithm developed by Lam 
and Goussis [14]. The algorithm examines the time 
scales of ODE systems and supports the separation of 
the biochemical network into fast and slow subsystems. 
This is achieved through the elimination of the detected 
quasi-stationary species and quasi-equilibrium reactions. 

The original CSP algorithm is implemented in the 
software COPASI [15] making it accessible to the scien- 
tific community. COPASI is a platform-independent, 
user-friendly software tool that allows easy access to 
powerful numerical methods for simulation and analysis 
of biochemical reaction networks. 

We apply the simplification strategy to two different 
systems to exemplify its use. Thus, as a simple system, 
we present the derivation of the Michaelis-Menten 
Kinetics. As a realistic case, we analyze the glycolysis in 
Saccharomyces cerevisiae [16] in three different dynamic 
regimes. We show that several variables can be elimi- 
nated still keeping the original dynamics intact. Further- 
more, regulatory mechanisms cause different players to 
participate with different relative importance in the 
dynamics of the system. 

Time Scale Separation Analysis 

In order to explain the basic notions of a time scale 
decomposition we start with a first-order kinetics sys- 
tem. Then, the differential equations describing the sys- 
tem dynamics y are linear: 



with constant and diagonalisable Jacobian J. By using 
the set of right eigenvectors A of J as the new basis we 
can decompose the Jacobian [17] and transform the ori- 
ginal equations into: 

x = A -1 -y, A = A 1 J A. 



The components of the transformed concentration 
vector x are called modes. Because A is a diagonal 
matrix of real or complex eigenvalues X t of J, the trans- 
formed ODE system is fully decoupled: 

dx l . . 

— = x l x l , i = 1, N. 
dt 

Thus, the modes x i ^ e ^t evolve independently of 
each other. The reciprocals of $R(A'): 

1 

have a dimension of time and are called time scales 
(TS). Ordering them w.r.t. magnitudes t x <t 2 < ... <t n 
leads to approximate speed ranking of the modes [14]. 
The modes corresponding to fast time scales (eigenva- 
lues with large negative real part) approach 0 very 
quickly and can be eliminated from the system for t » 
t m , where t m is a fast time scale. 

Two additional aspects are worth being emphasized 
here. First, although the transformed representation of 
the system dynamics in terms of modes provides a sys- 
tematic basis for the decomposition of the reaction net- 
work, it does not guarantee reducing the number of 
biochemical species or reactions in the system, since 
many different species might contribute to one and the 
same transformed equation. So, there is no straightfor- 
ward relation to reduction methods commonly used in 
biochemistry such as the quasi steady state (QSS) 
approximation or the quasi equilibrium (QE) 
assumption. 

An additional aspect of TS decomposition is that it is 
based on the local analysis of the system dynamics. For 
general nonlinear problems however the Jacobian is 
time-dependent. Its eigenvalues and eigenvectors change 
with time. Hence, in order to obtain a reasonable char- 
acterization of the systems dynamics the time scale 
decomposition has to be applied repeatedly at many 
time points through the evaluation time of the reaction 
system. 

2 Results 

2.1 CSP in COPASI 

Consider a system consisting of K biochemical reactions, 
the dynamics of which is determined by a system of N 
ordinary differential equations: 

^ = g(y(0) = E s ^ r (y) (D 

r=l 

here y is the A/-dimensional concentrations vector, s r (r 
= 1, . . . , R) are the Af-dimensional stoichiometric vec- 
tors and Fty) is the rate of the r-th reaction. 
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The main idea of the CSP method is to split the N- 
dimensional space of the vector g into two subspaces, a 
fast and a slow subspace: 



g _ gfast _j_ gslow 

In general, an A/-dimensional vector may be expressed 
in terms of any set of N linearly independent basis vec- 
tors (e.g. [17]). The objective of the CSP method is to 
express g in a new basis, one that is tuned to the 
dynamics of the system, where the fast and slow compo- 
nents evolve independently of each other. 

The subspace g fast relates to the fast time scales of the 
system. If its contribution is negligible (according to 
some error criteria), the original system (Eq. 1) simpli- 
fies to the system of the following differential algebraic 
equations (DAEs): 



dt 

fast ^ i 



i = M + l,...,N, 
1,...,M. 



(2) 



This DAE system does not contain the fastest time 
scales of the original model. Hence, it is much less stiff 
than the original system and can be simulated easily. 
Nevertheless such a simplification does not guarantee 
the reduction of the number of species and reactions 
("rear system reduction) as explained above. Therefore, 
a similarly important aspect of our strategy is the identi- 
fication of QSS metabolites and QE reactions by means 
of CSP data (see [14] and below). 
The Algorithm 

Let us differentiate equation (1) with respect to time 
and get the following form with the Jacobian J: 



3y 



(3) 



Now we (again) focus on the choice of an appropriate 
basis i = 1, . . . , N, but in contrast to the approach 
of time scale separation, these linearly independent vec- 
tors do not have to be eigenvectors of the Jacobian J 
(not even orthogonal). Then g always has the unique 
representation: 

N 

g = !>/'' 

i=l 

where a,/ is a so-called reaction mode, the amplitude 
/ is given by: 

R 

n Y ) = b i Qg = J2KF, i=l,2 N, 

r=l W 

B\ = b 1 O s r/ i = 1,2, . . . , N, r=l,2, R. 



The notation 0 abbreviates the standard scalar pro- 
duct, i.e. here b 1 O a ; = J2n=i tyitf- ^ ne set °^ ^ row 
vectors b l are the inverses of a r ; together they satisfy the 
following orthonormal condition: 

biOa^a], i,j= 1,2,...,N. 

The CSP provides an algorithm to determine the 
number of fast modes M, and to compute the sets of 
linearly independent a, and h\ such that 

M N—M 



z=i 



Differentiating equation (4) with respect to time we 
get: 



dt 



f=Bg, 



dB 



dt 



+ B J A, A = B" 



(5) 



(6) 



where A and B are matrices consisting of the column 
vector a, and row vector b l as basis vectors. 

For linear problems the matrix A is time independent. 
In this case the choice of eigenvectors of the Jacobian as 
the new basis leads to the diagonal matrix A (see 
above). The corresponding amplitudes f evolve indepen- 
dently of each other with their own characteristic time 
scale T/. For general nonlinear systems, however, A is 
time dependent and usually not diagonal. The CSP 
method provides an iterative procedure of refinement of 
basis vectors a/ and b ; . When recursively applied, the 
refinement procedure weakens the coupling between the 
M fast and the N - M slow amplitudes. The matrix A 
built from the final refined set of basis vectors is block- 
diagonal and the fast amplitudes are uncoupled from 
the slow ones approximately, so that the residual cou- 
pling can be neglected. 

The process starts with an arbitrary initial guess for 
the basis vectors a, and the assumption that the first M 
basis vectors span the M-dimensional fast subdomain. 
The corresponding time scales should be much faster 
than some characteristic time scale of interest, then M 
is to be selected to provide a gap between the slow 
(time of interest) and fast time scales: 



< e. 



(7) 



When for the final set of basis vectors the sum of M 
fast reaction modes falls below some user-specified 
threshold: 
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M 



£rel • Yj + £abs 



(8) 



these can be eliminated from the initial system (Eq. 1), 
because their contributions to g are negligible. As a con- 
sequence, the evolution of the reduced system depends 
on the slow modes only (see Eq. 2). 
Implementation 

The CSP algorithm was implemented as an integral part 
of the COPASI software in C++ and is freely available 
with the current releases of the package. In this imple- 
mentation, the CSP algorithm is applied to the models 
for which linear dependencies due to conservation rela- 
tionship are eliminated. This is achieved by the analysis 
of the stoichiometric matrix and is performed by 
COPASI automatically. 

The following CSP parameters have to be defined by 
the user. 

Intervals The user specifies the number of time points 
for which the CSP analysis is carried out by setting the 
time interval. The time interval should be large in com- 
parison with the user's time scale of interest. 
Ratio of time scale separation s This parameter speci- 
fies the gap between the time scales related to the fast 
and slow modes (Eq. 7). 

Error tolerance Absolute s^bs an d relative error e re \ are 
set to control when a fast mode is considered to be 
exhausted (Eq. 8). 

The CSP algorithm described above provides local 
information at certain time steps. To obtain global fea- 
tures of the system behavior the analysis must be per- 
formed at all points in the range of interest. For this 
purpose the CSP step involves numerical integration 
using the LSODA solver [18]. LSODA is part of the 
ODEPACK library [19]. It solves ODE systems with a 
dense or banded Jacobian when the problem is stiff, but 
it automatically selects between non-stiff (Adams) and 
stiff (BDF) methods. The Jacobian is generated 
numerically. 

When the CSP algorithm at time point t has been per- 
formed and both final refined sets of basis vectors a t -(£) 
and b\t) are available the M(t) is set to the number of 
fast exhausted modes and T M {t) is then the time scale of 
the slowest of fast reaction modes at time point t. 

The CSP output data (see below) can either be 
exported to a text file (save as Report in COPASI) for 
the use in other software (gnuplot, Octave etc.) or dis- 
played in the graphical user interface as tables. In this 
case a color coding is used where the numbers are addi- 
tionally visualized by different shades of color. This 
makes it easy to immediately spot e.g. the most impor- 
tant contributions to a specific mode for a large model 
(where the result tables are correspondingly large). 



We also use three dimensional bar graphs for visualiz- 
ing the matrices employing the qwtplot3D library 
(http://qwtplot3d.sourceforge.net) integrated in COPASI. 
These bar graphs can be turned and zoomed interac- 
tively. Furthermore single rows or columns of the 
matrix can be highlighted. An additional diagram shows 
the distribution of the time-scales of the different modes 
at chosen points of time (Figure 1). Applying the time 
slider in the graphical user interface it is very simple to 
switch between the results for different time points. 
Therefore the user can easily get an overview of the 
time-dependent changes of the time-scale separation. 
CSP Data used for model analysis 

The CSP algorithm supplies the modeler with local CSP 
output data [14] that relates the time scales to species 
and reactions of the original biochemical system. The 
data is computed by the help of the refined sets of basis 
vectors a,(£) and b\t). The user is provided with the 
CSP output at each defined time point during the inter- 
val of interest and can use it to reduce the model in a 
rational way. The CSP output data are displayed in 
COPASI in a number of matrices. Here we briefly 
explain the most important CSP output data which are 
available in COPASI: 

Time scales The analysis of time scales evolution can 
provide useful information about the system dynamics. 
The fast dissipative time scales relate to the eigenvalues 
of the Jacobian with large negative real parts. The explo- 
sive modes are associated with positive eigenvalues. 
Modes with equal time scales correspond to pairs of 
complex conjugate eigenvalues indicating oscillatory 
components in the system behavior. 



f 
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Figure 1 COPASI visualization of the time scale distribution. 

Full glycolysis model of Hynne et al. [16], [Glc x ] 0 = 14 mM, t = 25 
min. The coincident bars on the graph correspond to equal time 
scales. 
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Radical Pointer (RP) The CSP Radical pointer identi- 
fies the species for which the QSSA can be justified. 
Whenever the /-th diagonal element of m-th fast mode 
projection matrix Q m = a m b m is not a small number, 
species i is said to be a CSP radical 
Participation Index (PI) and QE reactions The relative 
level of participation of the r-th elementary reaction to 
the n-th CSP reaction mode can be represented by the 
mode participation index, P l r , defined as follows: 



BIF 



E \W\ + 

r=l 



where i - 1, . . . , N, r = 1, . . . , R. 
Importance Index (II) The relative importance of the 
contribution of the r-th elementary reaction to the rate 
of change of the z-th element of y can be represented by 
the importance index, l\\ 



R 

E 


s i,slow F| + 


Yerror 


r=l 







i = 1, . . . , N, r = 1, . . . , R. An effective stoichio- 
metric vector of the r-th elementary reaction, 
(I — Q(A4)) O s r is computed using the fast sub- 



r slow 



space projection matrix Q(M) = Qm - The reaction 

m=l 

with the largest l\ for the species y t is the rate control- 
ling reaction. 

CSP - based model reduction 

In this paragraph we summarize the most important 
steps in the reduction of the kinetic mechanism based 
on the results of the CSP algorithm described above. 
Model reduction is mainly the outcome of a sequence 
of QSSA for species and QEA for reactions which leads 
to the lumping or elimination of corresponding vari- 
ables. The QSSA identifies species whose production 
and destruction rate are in approximate balance. Mathe- 
matically it means that the right-hand side of the corre- 
sponding differential equation is zero. The QE 
assumption corresponds to reactions whose forward and 
reverse rates are nearly equal (see for instance, [20]). In 
either case an approximate algebraic relation (equation 
of state) is obtained between participating species. 

As described in the previous paragraph the CSP 
method provides the numerical data (RP, PI and II) that 
are an effective diagnostic tool allowing the detection of 
species which can be approximated by an equation of 
state, as well as the determination of the relative level of 
participation of distinct reactions to the modes. 



In contrast to the original CSP method [14] we intro- 
duce and use the "subspace" radical pointers and the 
"subspace" participation indices rather than the indivi- 
dual mode RP and PL This is based on the fact, that 
even though the matrix A (Eq. 6) built by the help of 
the final refined set of the basis vectors is block-diagonal 
and the fast modes are decoupled from the slow ones, 
the fast and slow modes could be coupled between 
themselves. So, it appears to be more reasonable to con- 
sider a projection of the CSP indices on the full fast and 
slow subspaces. 

We consider the sum of all CSP radicals as selected by 
M fast modes and define the species with the largest 
"fast subspace" radical pointers as QSS. Similarly the 
sum of Participation indices over all slow and fast 
modes should be considered separately in order to 
detect the fast reactions. The normed Pis over fast and 
slow subspaces are: 



PI 



fast 



M 
i=l 

1=1 



Pi! 



slow 



N 

E Pi 

i=M+l 

ep4 

1=1 



(9) 



We declare the reaction k as QE, if it is active in the 
fast and does not influence the slow space: 
pjfast ^ pjslow at all time po i nts> w here the CSP analysis 

was carried out. 

Practically, there exist only very few guidelines in the 
literature for deriving model simplifications based on 
the QEA and QSSA. Therefore, we would like to quickly 
summarize the procedure for the CSP-based model sim- 
plification: 

1. First, a time scale of interest should be selected. 
This can for instance correspond to the time resolu- 
tion of the experiment which is the basis for the 
model. The aim of the model simplification is to 
reduce all scales that are faster than this chosen 
scale. 

2. Second, user defined parameters have to be 
selected in COPASI as explained above. Since the 
CSP information will be available for every time 
interval and is the basis for the time-dependent 
model reduction, the time interval should be large 
enough in comparison with the time scale of 
interest. 

3. Third, performing the CSP and analyzing the 
results in order to find the QSS species and QE 
reactions. 

4. Fourth, solving the corresponding algebraic equa- 
tions of state and eliminating the respective variables 
from the reaction networks. The kinetic laws for 
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slow reactions should be modified by substitution 
with explicit expressions for CSP radicals. 
5. Fifth, parameter adjustment (e.g. by parameter 
estimation) with respect to quasi equilibrium con- 
stants of eliminated reactions in order to achieve the 
desired accuracy. 

It is worth mentioning that during the simplification 
the existing conservation laws have to be preserved. The 
algebraic equations should be solved under conditions 
that the equations of moieties are fulfilled. 

2.2 Application examples 

We have applied the method to two models. All infor- 
mations and scripts needed to reproduce the figures in 
this subsection are available in Additional Files 1, 2, 3. 
Michaelis-Menten Kinetics 

As in [10] we start our discussion with the simplest 
enzymatic reaction mechanism, the irreversible Michae- 
lis-Menten kinetics: 

Substrate Enzyme u Complex u Product Enzyme 

S + E <d> C -4 p + E. 

fe-i 

The model was build in COPASI and consists of the 
two reactions R x (5 +£ o C) and R2 (C -> P +E). In 
order to illustrate the handling of the CSP based model 
reduction we consider two limit situations for the 

k 2 

dimensionless parameters: S t = > 0 and 

fe-i 

M r = — > 0 (here E 0 and S 0 are initial enzyme and 

substrate concentrations, respectively). We used the fol- 
lowing CSP parameter values: s = 0.01, e rel = 10" 5 , and 
£ a bs = 10" 10 . In both cases, a clear time-scale separation 
occurs. 

(i) M r — » 0 (E 0 « S 0 ).* This is the standard situation 
for Michaelis-Menten kinetics. The motion on the fast 
time-scale is dominated by the complex C decoupled 
from the substrate S. On the slow part, the changes of 
substrate and complex are balanced. The quasi steady 
state assumption for complex C leads then to the 
Michaelis-Menten kinetic law. 

The CSP method allows the distinction between slow 
and fast modes (for times t >0.03). The Radical Pointer 
from the CSP data shows that the complex C dominates 
the fast mode. The contributions of both reactions to 
the slow and fast modes are comparable (see Figure 2, 
which displays the evolution in time of Radical Pointer 
and Participation Indices). Thus, the QSSA for the com- 
plex C is justified in this case. 

(ii) S t — > 0 (k 2 « k-i). This limit means that an equi- 
librium between the enzyme £, the substrate S and the 
enzyme-substrate complex C is established quickly. The 



slow step is the breakdown of C to produce the product 
P and the enzyme E. 

The CSP analysis leads to the occurrence of two inde- 
pendent dynamical modes. After the short transient 
phase (t < 0.006, when no reduction is possible) the con- 
tribution of C to the fast mode is larger than the one of 
S (nevertheless no real dominance occurs). Over the 
time the contribution of both variables becomes equal 
(Figure 2). Thus, the QSSA for complex C is incorrect. 

Nevertheless, there is a clear separation of reactions in 
the modes. The reaction R 2 of product formation domi- 
nates clearly the slow mode. Both reactions are active in 
the fast mode (see Figure 2). Thus, the reaction 

Ri : S + E <k- C is always practically in equilibrium and 

k-i 

the QEA for reaction Ri is correct and leads to a similar 
equations as for "standard" Michaelis-Menten kinetics 
(compare [21] and [10]): 

dP k 2 E 0 S , fe_i 

— = , with equilibrium constant K s = - — . 

dt K s + S H h 

The reader is referred to Additional File 1 for more 
details of the reduction of the Michaelis-Menten 
kinetics. 

Glycolysis in Saccharomyces cerevisiae 

We now use the CSP method to examine a more com- 
plex model for simplification purposes. We take the 
quantitative model of yeast glycolysis developed by [16] 
as an application example which has been also used 
before in similar studies [13,22]. 

The model is based on ODEs and consists of 24 reac- 
tions among 22 metabolites with a total of 59 kinetic 
parameters. The reaction scheme is depicted in Figure 
3. From the reaction stoichiometries two moiety conser- 
vations are derived: 

NAD + NADH = const; 
ATP + ADP + AMP = const. 

The complete model is available for download in 
SBML format at the BioModels database [23] (BIOMD 
61) or JWS online [24] (http://jjj.biochem.sun.ac.za/data- 
base/hynne/Hynne.xml), the latter version being used in 
this study. For model details the reader is also referred 
to [16]. However, there are some model properties we 
want to mention here explicitly. 

The model reproduces experimental data observed in 
intact yeast cells in a continuous-flow stirred tank reac- 
tor. Here, the mixed flow glucose concentration, [GlcJ 0 , 
is a bifurcation parameter which means that depending 
on its value the system behavior changes qualitatively. 
To be concrete, this glycolysis model exhibits two sta- 
tionary (<9.6 mM; 16.7 <[GlcJ 0 <18.5 mM) and two 
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(a) M r 0 (b) S t ^ 0 

Figure 2 Michaelis Menten model. Left: M r -> 0 (S 0 = 100; E 0 = 1 ; /c n = 1 ; = /c 2 = 100). Right: S t -> 0 (S 0 = 100; £ 0 = 100; ^ = = 100; k 2 = 
1). Time evolution of the Radical Pointer (RP) in the fast mode (top), Participation Indexes (PI) of reactions R ] and R 2 in the fast (middle) and slow 
mode (bottom). The RP of product P in the first case M r — > 0 is similar to RP of substrate S (the both lines overlaid). 



oscillatory state regimes (9.6 < [GlcJ 0 < 16.7; > 18.5 
mM). Please refer to Figure eight in [16] for the bifurca- 
tion diagram. Notably, the first oscillatory regime has 
not yet been observed in experiments. So, we consider 
this as an important model property. 
First step: CSP Analysis in COPASI When performing 
a model reduction analysis it is indispensable to deter- 
mine beforehand which properties of the system are to 
be maintained in the simplified model. We aimed at 
preserving (within an acceptable error range) the follow- 
ing features in order of priority: 

1. A Hopf bifurcation occurs at some value of [G\c x ] 

o 

2. Bifurcation points w.r.t. [GlcJ 0 change only little, 
i.e. different dynamic regimes (including the first 
oscillatory domain) appear at values of [GlcJ 0 close 
to the corresponding values in the full system. 

3. Steady state levels of metabolite concentrations. 

4. Periods of the oscillations. 

5. Amplitudes of the oscillations. 

We, therefore, perform the CSP analysis on the differ- 
ent dynamic regimes separately, i.e. using three different 
initial conditions for [GlcJ 0 , namely 9 mM (steady 
state), 14 mM and 24 mM (first and second oscillatory 
state, respectively). All other parameters of the model 
are taken as in [16]. 

For each CSP analysis we simulate the system for a 
time period from 0 to 100 min, thereby taking also the 
initial transients into consideration, and inspect 250 
time points along the trajectory which yields a time 
interval of 0.4 min. At each time point a full set of CSP 



data is computed. Example time course trajectories of 
the concentrations of ATP and NADH are shown in 
Figure 4. The CSP parameters Ratio of mode separation, 
Relative Error and Absolute Error are set to 0.99, le-3 
and le-4, respectively. 

In the following, we present the CSP output data 
(Time Scales, Radical Pointer, Participation Index, 
Importance Index and so on, see 2.1) one after the 
other. For each type of data, we point out the major dif- 
ferences between the three dynamic regimes which we 
interpret as glucose-dependent phenomena. If appropri- 
ate, special emphasis is given to time-dependent 
differences. 

Since the amount of data produced in this compre- 
hensive analysis exceeds the scope of the paper we pre- 
sent each CSP output data with compelling examples. 
The complete set of data is provided in Additional file 2. 
Time scales 

The full model exhibits in total 20 different time scales 
with values that span about seven orders of magnitude 
(from min to ms). Figure 1 shows the time scale distribu- 
tion (logarithmic values) of the full model exemplarily for 
[GlcJ 0 = 14 mM at time step 25. Notably, the time scale 
values change over time. In the steady state regime ([GlcJ 
0 = 9 mM), we observe two eigenvalue pairs corresponding 
to the 8 th and 9 th as well as 15^ and 16^ time scales that 
consist of complex conjugates (r 8 = r 9 , r 15 = r 16 ) indicating 
the system's intrinsic oscillatory vicinity. We see that the 
real part of these eigenvalue pairs become equal at a cer- 
tain point in time during the initial transient (Figure 5(a)). 

In both oscillatory state regimes, after the initial tran- 
sients, the values of time scales become oscillating and 
show in part substantial amplitudes which sometimes 
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also overlap with the values of adjacent time scales. As Number of fast modes M (Figure 5(d)): As explained 

an example, Figures 5(b) and 5(c) show the time evolu- above, each time scale corresponds either to a fast or 

tion of the 13 th to l& th time scales for [GlcJ 0 = 14 and slow so-called mode in the CSP analysis. Like the values 

24 mM, respectively. of the time scales the number of modes constituting the 



Surovtsova et al. BMC Systems Biology 2012, 6:14 
http://www.biomedcentral.eom/1 752-0509/6/1 4 



Page 9 of 16 




40 60 

Time [min 



(a) 9 mM 



(b) 14 mM 



97 98 99 

Time [min] 

(c) 24 mM 



Figure 4 Simulated time courses of [ATP] and [NADH] in the three different dynamic regimes at concentrations of [Glc x ] 0 from time f 

= 0 min to t = 100 min. In Figure (b), (c) the subinterval from f = 96 min is drawn to a larger scale. 



25 



2>20 
o 

0) 
| 

i- 5 



time scale 15 
time scale 16 
time scale 17 
time scale 18 



20 



40 60 

Time [min] 

(a) 9 mM 



80 



100 




96 97 98 99 

Time [min] 

(b) 14 mM 



100 



50 

45 

0)40 

05 35 
> 

0)30 

S25 
if) 

0)20 

10 
5 



30 
25 
20 
15 
10 

R 






0 20 40 60 80 100 

AAAAAA/ 




77 



96 



97 



98 

Time [min] 

(c) 24 mM 



99 



100 



Y 



Glcx6=9- 
GlcX0=14- 
GlcX0=24- 




(d) Fast modes 



Figure 5 Time evolution of the time scales 15 to 18 in the three dynamic regimes (a)-(c) at concentrations of [Glc x ] 0 from time t = 0 
min to t = 100 min. In Figure (b), (c) the subinterval from f = 96 min is drawn to a larger scale, (d) Time evolution of the number of fast 
modes M in the three different dynamic regimes. 



Surovtsova et al. BMC Systems Biology 2012, 6:14 
http://www.biomedcentral.eom/1 752-0509/6/1 4 



Page 1 0 of 1 6 



entire fast or slow subspace changes over time. Since for 
model reduction only the fast modes are relevant we 
focus on these. Initially, all three dynamic regimes show 
seven fast modes. In the steady state regime, after a 
highly variable transient, M settles to 17. In contrast, M 
varies between 7 and 9 for the first and between 9 and 
10 for the second oscillatory regime. Consequentially, 
we do not fix the number of fast modes in our CSP ana- 
lysis but rather take their varying number over time into 
account in search for QSS metabolites (see RP) and QE 
reactions (see PI). 
CSP Radical Pointer 

Figure 6 shows how Radical Pointers are visualized in 
COPASI. Five metabolites (BPG, GAP, PEP, F6P and 
NAD) are fast in all of the three dynamic regimes. 
CSP Participation Index (PI) 

When comparing the normed sum of Pis for the three 
different regimes, four different categories of reactions 
can be identified depending on their respective Pis, e.g. 
the reaction can always be classified as fast or it changes 
its role between regimes. A heuristic threshold value 
based on our analysis and experience is chosen. Thus, if 



the normed sum of Pis over all fast modes exceeds 0.7, 
the reaction is defined as fast. 

1. "fast -fast -fast": vGAPDH, vlpPEP, vPK, vPGI, 
vALD, vTIM and vAK are fast in all regimes. These 
reactions, therefore, may be approximated as QE and 
eliminated in a simplified model. Not surprisingly, the 
known fast reactions vPGI and vTIM turn up in this 
group. Interestingly, the group also contains all reactions 
that either produce energy or redox equivalents, i.e. 
ATP and NADH, respectively. Obviously, especially in 
case of reactions being at the edge of the threshold, 
model reduction still has to be done with care. 

2. fast - slow - slow": vHK, vPFK, vPDC, glycerol pro- 
duction, glycogen production, and ATP consumption 
are reactions that belong to this group which switch 
from fast to slow with increasing [GlcJ 0 . These reac- 
tions (except vPDC) share the property of consuming 
energy and redox equivalents, i.e. ATP and NADH, 
respectively. The continuous flow transport reactions 
between the outside and the chemostat (vinCN, vinGlc) 
as well as vlacto also belong to this group. 
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3. "fast - slow - fast": vADH and the transport reac- 
tions across the cell membrane (vGlcTrans, vdifACA, 
vdifEtOH, vdifGlyc) behave differently from all others as 
they are fast for low and high concentrations of [GlcJ 0 . 
Participation in slow modes seems to be limited to the 
first oscillatory regime. 

4. "slow - slow - slow"-. All reactions from the chemo- 
stat to the outside (voutEtOH, voutGlyc, voutACA) are 
slow in all regimes. 

A typical example of time evolution of the normed Pis 
for each class of reactions is given in Figure 7. 
CSP Importance Index (II) 

The majority of reactions exhibit significant importance 
on a number of metabolites (Normed Importance Index 
>0.1). Exceptions are vPGI, vALD, vTIM, vlpPEP, vPK, 
vconsum, vAK and vdifACA, where Importance Indices 



are of values less than 0.1 for all metabolites. The weak 
importance of the first five reactions (already indicated 
as QE by the normed Pis) further confirms that they 
may be removed from the model. In some cases, the 
importance index changes in between regimes, depend- 
ing on [GIcJq. Examples for glucose-sensitive impor- 
tance are vinGlc (important at low and unimportant at 
high glucose concentrations), vHK, vPFK and vGAPDH 
(unimportant at low and important at high glucose con- 
centration). Obviously, the importance index gives simi- 
lar information as control coefficients derived from 
MCA, a fact that we studied and verified (data not 
shown). However, the CSP lis give a richer picture of 
the control distribution compared to MCA. 
Second Step: Model Reduction Based on the time scale 
separation analysis we suggest four steps to derive a 
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simplified minimal model. A short description is given 
in the following. For any detail the reader is referred to 
Additional file 3. Each simplification step concerns a 
subset of the original model scheme which we call Mod- 
ule, hereafter. 

Module 1. QEA for vPGI, AE for F6P. The normed 
PI revealed that PGI can be approximated as QE and 
the Radical Pointer of the 5-th fast mode identifies F6P 
as CSP Radical, for which the algebraic equation holds 



^PGI 



F6P 
G6P 



So, in order to eliminate F6P from the system and to 
lump PGI together with PFK we need to modify the 
chemical equation of the PFK reaction to 

G6P + ATP -> FBP + ADP 

as well as the kinetic rate law to 



V 5m • (K PGI • G6P) 2 



K 5 



/ATP 

1+K 5 - 

\AMP 



+ {K PGI • G6P) : 



Module 2. QEA for vALD and vTIM. The normed PI 
revealed that vALD and vTIM can be approximated as 
QE, for which the equations hold 



GAP • DHAP 



FBP 



and Ktim 



GAP 
DHAP' 



The metabolites which are either substrate or product 
of the two reactions are FBP, DHAP and GAP. The lat- 
ter is identified as CSP Radical (see Radical Pointer of 
the 2-nd fast mode). In order to lump vALD and vTIM 
together we introduce a pool metabolite which we name 

trioseP = GAP + DHAP + FBP 

and express any of the three metabolites in terms of 
trioseP. The new chemical equations of the associated 
reactions are: 



PFK 
GAPDH 
Glycerol branch 



G6P + ATP -> 2 • trioseP + ADP, 
trioseP + NAD -> BPG + NADH, 
trioseP + NADH -> Glyc + NAD. 



Module 3: QEA for vlpPEP. The equilibria for the 
vlpPEP reaction is expressed as: 



^PEP 



BPG • ADP 
PEP • ATP ' 



BPG is identified as CSP Radical in the first mode and 
at the same time PEP in the 4-th mode. Again, we intro- 
duce a pool metabolite 



BPG-PEP = BPG + PEP 

and reduce the vlPEP reaction from the network. The 
new chemical equations of the associated reactions are: 

GAPDH : trioseP + NAD -> BPG_PEP + NADH, 
PK : BPG_PEP + 2ADP -> Pyr + 2ATP. 

Module 4: QEA for vPK. The vPK reaction is mod- 
eled as irreversible. So, the QEA leads to its lumping 
together with vPDC and to eliminating pyruvate from 
the network. The new chemical equation for vPDC is: 

BPG-PEP + 2 ADP -> ACA + 2 ATP. 

In summary, after these four simplification steps the 
full model (original values in parentheses) has been 
reduced eventually to 17 (22) species and 19 (24) reac- 
tions with a total of 43 (59) parameters (the reduced 
reaction network is depicted on the Figure 8). 
Third step: Parameter adjustment and verification of 
the reduced model Due to the fact that the meaning of 
parameters has been changed in the course of model 
reduction these parameters (e.g. K4eq) need to be 
adjusted in order to obtain the full original behavior. 
This can be simply achieved by parameter scanning 
around the initial value. It is worth emphasizing here, 
that not all parameters have to be refitted, only the ones 
that result from the simplification of the lumping terms 
(e.g. quasi equilibrium constants resulting from the 
QEAs). 

Finally we evaluate the reduced model by comparing 
its dynamic properties with the ones of the original full 
model. Comparative simulations are shown in Figure 9 
and reveal that the reduced model captures the essential 
dynamics of the full model quantitatively very well - 
except for the amplitudes and the exact location of the 
bifurcation points for the first oscillatory regime. This 
discrepancy is of (only) quantitative nature and it does 
not occur if the full model is reduced by just three reac- 
tions (instead of five) as presented in Additional file 3. 

3 Discussion and Conclusions 

In this paper, we have presented a strategy for model 
simplification and reduction based on the CSP method. 
For this purpose and in order to make the method pub- 
licly available we implemented the original CSP algo- 
rithm in the COPASI software. 

The CSP method is restricted to ODE models. Pre- 
viously described simplification routines based on CSP 
mainly focus on the conversion of ODEs into DAE sys- 
tems. In contrast, we use the CSP method to simplify 
models by lumping those reactions together that could 
be identified as being in QE. In addition, algebraic equa- 
tions are used for species that are identified by Radical 
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Figure 8 Modified part of the reaction scheme for the reduced glycolysis model of S. cerevisiae. 



Pointers. Accordingly, we redefine chemical equations 
and kinetic rate laws of affected reactions. We demon- 
strated the usability of this approach using the COPASI 
implementation of the CSP method for a simple one- 
enzyme reaction and for a rather complex model of 
yeast glycolysis [16]. 

The time scale separation analysis of the glycolysis 
model revealed five reactions (vPGI, vALD, vTIM, 
vlpPEP, and vPK) for which the simplification strategy 
can be applied. We demonstrated that the resulting 
reduced model is capable of maintaining characteristics 
of the full model within an acceptable error range: 

(i) same dynamic regimes, e.g. Hopf bifurcation point 
at [Glc Jo = 18.5 mM; (ii) similar steady state levels of 
metabolite concentrations; (iii) similar periods for both 
and amplitudes for the second oscillatory regimes. 

Studying different dynamics underlines again (as in 
[11]) the importance of time-resolved analyses since the 
contribution of the players in the system may vary over 
time and in between different dynamical regimes. This 



is ignored if either steady state data (or single time 
point data in general) or single dynamic regimes are 
studied. 

Compared to our previous work on the ILDM method 
[10,11] - or the ILDM method in general - the CSP 
allows a more straightforward interpretation of its 
results with respect to the identification of QSS species 
and especially QE reactions. In addition, the Importance 
Index of CSP allows to analyze the impact of individual 
reactions on the dynamics of the species in the system. 

An interesting outcome of our analysis is that it is 
possible to follow the general inherent temporal organi- 
zation of the entire system when analyzing the distinc- 
tive time scales. Thus, we could observe that for the 
second oscillatory regime, all time scales oscillate in 
phase, partially overlapping each other which indicates 
that the whole system shows slower or faster dynamics 
in the course of a period. 

Moreover, the number of fast modes changes over 
time and is also different for different dynamic regimes. 
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Figure 9 Time courses of [ATP] and [NADH] in the two oscillating regimes at concentrations of [Glc x ] 0 as indicated. The upper diagrams 
show the simulation of the full model, the lower the ones of the reduced system. 



Both factors prohibit the use of a fixed number of 
modes for time scale decomposition. 

Furthermore, we suggest that the results of the CSP 
analysis can also be used for studying the relative 
importance of different reactions for the dynamics of 
the system. As an example, we observed that the overall 
participation of PFK in the slow modes increases with 
increasing glucose levels. In a simple way, this may be 
explained by the increasing energy charge (ATP concen- 
tration) which inhibits the PFK. Therefore, the relative 
importance of the PFK to the slower modes of the sys- 
tem increases. 

Another beneficial result of the simplification process 
is of course that the number of system parameters is 
considerably reduced, especially concerning parameters 
which are involved in processes on a faster time scale 
than the time scale of interest which are then usually 



hard to identify. Therefore, using this process less sys- 
tem parameters will be unidentifiable. 

Our study is not the first trying to reduce the original 
glycolysis model by [16]. [13] analyzed exclusively the 
limit cycle of the second oscillatory regime ([GlcJ 0 = 24 
mM) employing CSP without taking into account transi- 
ent behavior. In contrast, we analyzed the model with 
original initial values taking into consideration also the 
initial transient time period. In addition, there are major 
methodological differences. First, our approach focuses 
on simplifying the underlying biochemical reaction net- 
work rather than on approximating the ODE system 
with a DAE system. Second, we do not fix the number 
of fast modes. Third, we compute the normed sum of 
Pis over the entire fast subspace in order to justify QEA. 

A completely different approach was taken by [22]. 
Their sole criterium for the reduction was the 
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fulfillment of a Stuart-Landau equation which is in prin- 
ciple only valid in the vicinity of a Hopf bifurcation and 
therefore does not offer a general strategy for system 
reduction. 

Obviously, there are some relations between CSP out- 
put data and sensitivity analyses like metabolic control 
analysis (MCA). Learning e.g. about the impact of indi- 
vidual reactions on systems properties like dynamics 
could in principle also result from sensitivity analyses. 
We did a preliminary comparison of the results of our 
CSP analysis and a conventional MCA for the steady 
state. This resulted in a similar global picture, but the 
CSP gave a more fine-grained picture w.r.t. the relative 
importance of reactions on species. In addition, the 
time-resolved analysis for oscillations is not possible 
with MCA. 

With all the mentioned benefits of using CSP for sys- 
tems analysis, there are also problems and limitations 
arising from this approach. We employed several heuris- 
tic thresholds for the discrimination of the reactions and 
species mainly contributing to the fast subspace of the 
system. These were based on our experience and 
obviously, this might not be optimal for arbitrary sys- 
tems. Thus, other systems might demand slightly altered 
thresholds. This is underlined by the fact that we 
observed one reaction - AK - that in principle fulfilled 
all of our criteria for elimination, but in the end, it 
turned out to be impossible to eliminate from the sys- 
tem without introducing a large error. Therefore, it is 
always important to carefully check the behavior of the 
reduced system. The CSP can only support this process 
in a rational way, but does not allow for a fully auto- 
mated analysis. 

Even though, accordingly, scientists will always have to 
be on top of this method, it would be useful to support 
the reduction of the system in a stronger way than just 
providing the CSP. A semi-automated reduction which 
then quickly allows to be checked for error compared to 
the original model would reduce workload considerably 
and is currently planned to be included in the software. 
An additional planned extension of the software is the 
support of different compartment sizes (if multi-com- 
partment models are analyzed) which is currently not 
the case. 

All in all, we were surprised that taking into account 
different dynamic regimes only allowed the elimination 
of 5 reactions and 5 species of the glycolysis model 
which is considerably less than previous attempts that 
focused on particular regimes. This once again supports 
the view that it is crucial to define which systems beha- 
viors should be reproduced by the simplified model 
before entering reduction strategies and these initial 
decisions might result in different models in the end. 



Additional material 



Additional file 1: Michaelis Menten Kinetics: Includes a more detailed 
simplification procedure of Michaelis Menten kinetics. 

Additional file 2: CSP output data for glycolysis model: Includes the 
complete set of CSP output data (time resolved TS, RP, PI and II) for the 
glycolysis model. 

Additional file 3: simplification of glycolysis model. Additional 
material related to the simplification of the glycolysis model. This 
includes a list of original and modified reactions, kinetics laws and 
parameters. 
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